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Abstract. Smoothing is omnipresent in astronomy, because almost always measurements performed at discrete 
positions in the sky need to be interpolated into a smooth map for subsequent analysis. Still, the statistical 
properties of different interpolation techniques are very poorly known. In this paper, we consider the general 
problem of interpolating discrete data whose location measurements are distributed on the sky according to a 
known density distribution (with or without clustering). We derive expressions for the expectation value and for 
the covariance of the smoothed map for many interpolation techniques, and obtain a general method that can be 
used to obtain these quantities for any linear smoothing. Moreover, we show that few basic properties of smoothing 
procedures have important consequences on the statistical properties of the smoothed map. Our analysis allows 
one to obtain the statistical properties of an arbitrary interpolation procedure, and thus to optimally choose the 
technique that is most suitable for one's needs. 
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1. Introduction 

A common problem in Astronomy is the smoothing of ir- 
regularly sampled data. In general, this happens when one 
can make measurements on discrete points of a quantity 
that has some astrophysical relevance. In most cases, the 
"quantity" is a continuous field, i.e. a smooth function on 
the sky. In this case it is reasonable to try to reconstruct 
the field by interpolating the discrete measurements ob- 
tained. 

Suppose, for example, that we are interested in mea- 
suring the column density of a nearby molecular cloud. 
A good estimate of the cloud column density can be ob- 
tained, for example, using the infrared color of background 
stars observed through the cloud (see, e.g., Lombardi & 



Alves 2001). This way, we can obtain reliable measure- 
ments of the "column density field" at the discrete points 
corresponding to the star positions. Finally, we interpo- 
late the various measurements and obtain a smooth map. 
Similar situations are often encountered in different fields 
(e.g., weak gravitational lensing, peculiar velocity field of 
galaxies) . 

Interpolation and smoothing are ubiquitous in 
Astronomy, and indeed many papers have been devoted 
to the study of the effects of interpolation in particular 



the study of interpolation are possible, depending on the 
type of problem considered. In this paper, in particular, 
we will study the statistical properties of several interpola- 
tion techniques. Since the locations on the sky where mea- 
surements are performed cannot generally be chosen in 
Astronomy (cf. the example of reddening of stars above), 
we will carry out an ensemble average over the measure- 
ment locations. More precisely, we will assume that the 
measurements are randomly distributed on the sky fol- 
lowing a known scheme and we will carry out a statistical 
analysis on this sample. Note that there is a complete 
freedom on the spatial distribution of measurements on 
the sky, that can be correlated (for example, clustered) 
or can follow a non-uniform density p{x). The ensemble 
average, which has already been carried out under sim- 
plifying hypotheses of uncorrelated points and uniform 



density in earlier papers (e.g.. 


Lombardi & Bertin 


1998; 


van Waerbeke 


2000|; Lombardi & Schneider 


2001, 


here- 


after Paper I; Lombardi et al 


2001; 


Lombardi & Schneider 



analyses (see, e.g. Rybicki & Press 1992). However, it is 
important to observe that many different approaches to 
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that generalize the particular configuration considered. 

In previous works (see in particular Paper I and II) we 
have focused on a widely used smoothing method. Here, 
in contrast, we will keep the discussion much more gen- 
eral and consider some wide classes of interpolating tech- 
niques. This alternative approach has a number of advan- 
tages: (i) The results obtained are very general and can be 
applied to several smoothing techniques ( "proof reusabil- 
ity"); (ii) The analytic discussion is kept at a very simple 
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level; (iii) The results obtained are more general, since al- 
low for a non-uniform density of locations and correlation 
on the locations; (iv) The properties of a smoothing tech- 
nique can be predicted in advance without the need of long 
calculations. This last point, in our opinion, is particularly 
important for practical applications. In fact, with the aid 
of the results obtained in this paper, several smoothing 
techniques can be easily compared, which allows one to 
choose the interpolation method more suitable. We stress, 
however, that some specific results cannot be derived using 
the techniques described in this article and need a more 
specific analysis as done in Paper I and II. 

The paper is organized in two main parts. In the first 
part. Sect. ||, we classify the smoothing methods and we 
show the general results associated to each interpolation 
family. In the second part. Sect. ^, we illustrate the results 
obtained earlier by considering several common smoothing 
techniques. Finally, in Sect. ^, we briefly draw the conclu- 
sions of this analysis. 

1.1. Notation 

Let us suppose to have a set of pairs {(xi,?/^)}, where 
Xi G X, called "locations," belong to a real vector space 
X (typically, X = R", with n = 1, 2, or 3), and y, G Y, 
called "values," belong to a field Y (in this paper, we will 
assume for simplicity K = R). The spatial interpolation 
problem consists in finding a way to obtain approximated 
values for a generic x G X. Formally, an interpolation 
procedure is a function S that maps a point x € X and 
an unordered set of couples {xi, yi) G X xY into a value 
y^Y: 



y = S{x■,{{x^,y^)]) . 



(1) 



Note that the number N of couples {xi, yi) is not fixed a 
priori (in particular N could vanish). In the following we 
will use as synonymous interpolation procedure, smooth- 
ing technique, interpolator. 

As an example, let us consider a simple interpolator 
where the smoothed value at x is obtained as a weighted 
sum of the values {yi}. More precisely, we define S as 



S{x; {{xi.yi)]) = 



(2) 



mean and covariance matrix (taking {ei} as a multidimen- 
sional random variable) proportional to the identity: 



-0, 



(4) 



Note that the second relation of Eq. (Q) states the so- 
called "statistical orthogonality" of the measurement er- 
rors; this property is trivially satisfied if {e^} are indepen- 
dent random variables with fixed variance (this is a 
good approximation for many astronomical observations). 
The subject of our study will be the average value {f{x)) 
of f{x), where f{x) is the interpolated value of / at x: 

f{x) = S{x; {{x^,y^)}) = S{x;{{xiJ{xi) -l-e,)}) . (5) 

We will also investigate the covariance (or two-point cor- 
relation) of the map f{x), defined as 

Cov(x,x';/) = ([/» - (/»)] [fix') - (/>'))] 

= (/(x)/V)) -(/»)(/>'))• (6) 

A word of explanation is needed regarding averages. 
Averages of f{x) and f{x)f{x') are carried out both with 
respect to {ci} and to {xi}. As we have anticipated above, 
averages with respect to {x^} are carried out assuming 
that the locations are random variables distributed with 
density p{x) over the set X. For example, if we assume 
that the locations are uncorrelated and that the density 
is constant over the field (this is usually referred to as a 
homogeneous Poisson process or as compete spatial ran- 
domness), then the number N of object locations inside a 
field A (Z X finite area /^(A) follows a Poisson distribu- 
tion with average p^,{A): 



Pn{N) 



[p^^{A)] 

Nl 



N 



(7) 



Under the same hypotheses, the N locations {xi} are, 
then, uniformly distributed inside A. In more general cases 
(in particular, in case of correlation), it can be non-trivial 
to specify exactly the probability distribution of points. 
However, fortunately for the following discussion we need 
only the density p{x) and an algorithm to randomly gen- 
erate the locations (as discussed by Scott ct aL| 1954; see 



where a is a fixed real number, 'i'his interpolation tcch- 
nique will be often used in this paper to illustrate with an 
example some general, abstract results. 

In this paper we will study some statistical properties 
of various interpolation procedures assuming that the lo- 
cations {xi } are random variables distributed with spatial 
density p{x) (with or without clustering), and that the 
values {yi} are associated to the locations through the 
relation 



also, e.g., the hierarchical model described by 
Peebles|[l97^). 



50neira & 



= fi^i) + ei 



(3) 



Here / : X — > y is a known function and {e;}, representing 
measurement errors, are random variables with vanishing 



The probability distribution for the other random vari- 
ables considered here, namely {e^}, need not to be fully 
specified. Actually, for our purposes, we just need to spec- 
ify that these variables, representing measurement errors, 
have vanishing mean, fixed variance, are orthogonal [i.e., 
satisfy Eq. (^], and are independent of the locations {xi}. 
In the following calculations, we will carry out first the av- 
erage over the measurement errors {ci}, and then the one 
over the locations {xi}. 

We finally note that in this study we allow for cases 
where the smoothing procedure S cannot be defined. For 
example, if the density p is small, we might end up with 
a configuration without locations Xi close enough to x. 
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Some interpolation procedures then might not be applied. 
In the following, we adopt the convention of discarding 
in the ensemble average the configurations {xi } for which 
the interpolation procedure (|l|) is not locally defined. For 
example, when evaluating the average (/(a;)), we discard 
configurations that make S not defined at x; similarly, for 
expressions such as we discard configurations 

for which S is not defined at x or at x' . Note that we 
assume that the applicability of the smoothing technique 
depends only on the locations {xi} and not on the values 
{Ui}. We call Po{x) the probabiHty of having a configura- 
tion {xi} for which the estimator is not defined at x. For 



example, some smoothing techniques (see Sect. 2.6) are 
not defined if there are no locations inside a given subset 
tt; in this case we find (in case of vanishing correlation, i.e. 
for a Poisson distribution of locations) 



Po{x) = exp 



p{x) dx 



(8) 



Similarly, we define Po{x,x') as the probability that the 
smoothing is not defined either at x or at x' (or at both 
points). If, again, the smoothing technique is defined at x 
only if there is at least a location inside a given subset tt, 
and is defined at x' if there is a location inside the subset 
tt', then 



exp 



it\tt' 



p{x) dx 



exp 



7r'\7r 



p{x) dx 



exp 



(7r\ir')U(ir'\7r) 



p{x) dx 



(9) 



Note that Po{x,x) = Po{x) and that Pq{x)Po[x') < 
Po{x,x') < Po{x) + Po{x'). If an estimator is always de- 
fined, as for the interpolator (||), we just set Po{x) = 
Poix,x') = 0. 



2. Classification of interpolation procedures 

Very little could be said at this stage about the average 
of f{x) or its covariance because the problem in the for- 
mulation of Sect. 



1.1 



is by far too general. However, as we 
will see below, additional hypotheses allow us to obtain 
various results. Some of these hypotheses are quite obvi- 
ous, in the sense that we could hardly imagine a sensible 
interpolation procedure that does not satisfy them. For 
example, if all values are the same, say yi — 1, we expect 
to have always y = 1 in Eq. (|l|) (this property is discussed 
in Sect. 2.2). Still this "obvious" hypothesis allows us to 



obtain some non-trivial results. 

In this section we will consider some criteria used to 
classify interpolation procedures, with the aim of intro- 
ducing a standard terminology and, at the same time, of 
deriving a number of useful results. 



2.1. Linearity 

Some of the most interesting interpolation procedures are 
linear functions of the data values {yi}. We stress that the 
linearity is on the values and not on the locations {xi}. 
Formally, the linearity is expressed by the relation 

S{x; {{x^,ay^ + /3y-)}) = aS{x; {{x^,y^)}) 

+ PS{x;{{x,,yl)}) , (10) 

where a and (3 are arbitrary real numbers. Given a loca- 
tion configuration {xi}, we define the j-th weigJit of S as 



Wj{x]{xi}) ^ S{x\{{xi,5ij)]) 



(11) 



Using Eq. (|ll|) , we can write any linear interpolation pro- 
cedure in the form 



N 



S{x;{{x^,yi)}) ^^Wj{x]{xi})yj 



(12) 



For example, the smoothing technique (g) is manifestly 
a linear function of the values {yi}, and thus is a linear 
interpolator. Actually, it is already written in the form 
(ph with 



Wj{x;{x,}) = — — 



(13) 



Since the discussion of linear smoothing techniques is 
quite lengthy, we split it into two subsections, correspond- 
ing to the study of the bias and the covariance of the in- 
terpolated map. 

2.1.1. Bias 

Linearity is very convenient for a study of the bias of the 
smoothing, since it allows us to write the average value of 
a smoothed function / as 



(/(a;)) - j^f{x)K{^\^)p{^)^i 



(14) 



In other words, the bias of the estimator is completely 
described by the kernel function K{x\x). Note also that 
the only statistical property of measurement errors {e^} 
that have been used to derive Eq. ( p^ is (e;) = 0; again, 
this is due to linearity. Almost all interpolators considered 
in this paper are linear. 

In order to show Eq. (|lj), we consider the average 
{f{x)'j and use the linearity relation (10): 



(fix)) = (S{x;{{x,,f{x,) + ei)}) 



N 



= {J2[.f{x,) + e,]S{x;{{x,,S,,)}) 



(15) 



Here Sij is the Kronecker symbol. We consider now the 
average on errors, and note that the ej can be dropped 
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from the previous equation because {tj) = 0. Let us now 
decompose / in the form 



f{x) = / 5{x- x)f{x) dx 
Jx 



(16) 



where d is Dirac's distribution. Inserting this expression 
into Eq. (|l^) we obtain 

(/») = Q^dxfix)J2Six,^x)S{x-,{{x.,S,,)})^ 

= f f{x)K{x;x)p{x)dx , (17) 
where K{x] x) is given by 
1 / ^ 

K{x]x) = -^{YlHxj - x)S{x;{{x^,5ij)}) 



p{x) 



S{x; {{x^,S{xi - x)}) 



(18) 



Note that we have defined K{x; x) with a factor l/p{x) for 
further convenience. Equation ( p7[ ) is precisely Eq. (p^. 

The practical evaluation of the kernel K is very im- 
portant for a statistical study of the interpolation tech- 
nique. Indeed, this kernel controls the relationship be- 
tween the expected interpolated map and the original, un- 
known field, and thus is useful, for example, to perform 
a comparison between the observations and some theo- 
retical expectation. In order to evaluate K, we first ob- 
serve that almost all linear interpolation procedures have 
a simple property: The interpolated value at x does not 
strongly depend on values that are far away from it. In 
this case, the kernel K associated with a linear interpola- 
tor can be obtained using a numerical technique that we 
now describe (a proof will follow). The procedure used in 
order to determine K{x;x) is described in the following 
items: 

1. Choose a large subset A C X that contains both x and 

X. 

2. Generate N, the number of objects inside A, according 
to the scheme chosen for the distribution of the loca- 
tions (for example, if the locations are independent and 
uniformly distributed with density p, then N follows 
the Poisson distribution given in Eq. (Q)). 

3. Generate the N locations {xi} inside A following the 
scheme chosen (for example, uniformly if the density 
is constant and there is no correlation). 

4. Assign a vanishing value = to each location Xi. 
Assign a value 1 to an extra location at x. 

5. The configuration considered is composed by the ob- 
ject at X and the N objects at {xi}; hence the config- 
uration is T = {(a;, 1)} U {{xi, 0)}. 

6. Evaluate S{x] T) if this function is defined at x] oth- 
erwise arbitrarily assign a vanishing value to the inter- 
polated value at x. 



7. Generate several configurations by repeating points 2- 
6, and evaluate the average of ^(x; T). This average, 
multiplied by [l — ^0(2;)] , is an estimate of K{x; x). 

This procedure will be used in this paper not only to de- 
rive numerically the kernels of several interpolating pro- 
cedures, but also to obtain some general properties of in- 
terpolators. 

A proof of this practical method can be obtained as 
follows. From Eq. (|l^) we see that, in principle, K could 
be evaluated by performing an integration over the prob- 
ability distribution function for {xi} using the whole set 
X. In practice, this method could never be applied in nu- 
merical studies because of the presence of S Dirac's dis- 
tribution and because of the large number of locations to 
be generated (actually, since X is not bounded we should 
generate an infinite number of objects). However, a dif- 
ferent approach is feasible. First, we use a Monte-Carlo 
integration, that is, we generate a set of locations {xi} 
according to the expected probability distribution. If we 
assume that the smoothing is weakly dependent on values 
Hi whose locations Xi are far away from x, we can safely 
generate points inside a subset A d X that abundantly 
contains x and x. Then, in order to solve the problem 
with Dirac's delta, we approximate this distribution with 
a top-hat function 11(5; — x^) centered on x and with area 
s. In other words, we take II(x — Xi) to be 1/s if Xi falls 
inside a region of area s around x, and we take other- 
wise; we will eventually let s go to zero. Since s is taken 
to be small, in most cases all locations {xi} will not be 
close enough to x and thus all values {yi} will vanish. In 
such situations, since the interpolator is linear, we obtain 
a vanishing value at x. In a few cases, however, we ex- 
pect to have configurations with a single point inside s, 
and thus with value 1/s; such configurations have proba- 
bility p{x)s as s — > 0. When we now take the average of 
S, cases where all values vanish do not contribute to the 
average because the interpolated value vanishes as well; 
other cases, instead, contribute with a term proportional 
to 1/s (the value corresponding to the location around x). 
Equivalently, we can estimate the relevant average forcing 
a point inside the top-hat close to x and multiplying the 
resulting average by the probability that this happens, i.e. 
p{x)s. Since the interpolator is linear, this is equivalent to 
forcing a point with value 1 at x, taking the average, and 
multiplying the result by p{x)] this last factor, however, 
actually cancels with the term l/p{x) used in the defini- 
tion of K{x; x) [cf. Eq. ([l^)]. This proves the correctness of 
the procedure to obtain K described in the points above. 

In our discussion we have implicitly assumed that the 
estimator S is always defined, i.e. that Pi^{x) — 0. If this 
is not the case, we must slightly modify the method de- 
scribed above. In particular, we need to distinguish be- 
tween four probabilities: Ps,i, the probability of having 
a point inside s and the smoothing procedure defined 
at X, Py.i: the probability of having no point inside s 
and the smoothing procedure defined at x, and the two 
symmetrical probabilities Psfi and P^fi for configurations 
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where the interpolator is not defined at x. Note that Using this expression in the product (^f{x)f(x')'), we ob- 
-Ps.o + Pfifi = Po{x) as defined in Sect. LI. We then find tain 
that the procedure described above can be still applied if 
Po{x) 7^ 0, provided that two changes are made: 



(/»/>')) = (E[/(^^-) + ^^] [f(^f) + 'r] 
j-j' 

X S{x; {{xi,Sij)})S(x'; {{xi' ,Sej')}) 



— We need to discard, in the ensemble average of points, 
configurations for which S cannot be evaluated at x. 

— We need to multiply the average of S obtained from 
Monte-Carlo integration by the probability that, in the , 

ensemble average, a point falls inside s. In the gen- ~ \z_^[fi^j)fi^j)^'^ ]S{x; {{xi,Sij)}) 

eral case Po{x) ^ considered here, this probability is 
given by Ps,i/[s{Ps,i + P/,i)] . 



X S{x'; {{x,>,5i'j)}) 

Both problems can be solved at the same time using the 

following procedure. First, we note that the if we set y = _l_ / f[xj)f{xj')S{x- {{xi %)}) 

for cases where S cannot be evaluated at x, we are basi- S^i' 
cally multiplying the average by a factor Ps,i/{Ps,o+Ps,i)- \ 
On the other hand, we know from the previous analysis x «S'(a:'; {(a^i', <^i'j')}) / (23) 



that Pe,o + PsA ~ p{x)s. Hence, we can recover the correct 

factor by setting w = if S" is not defined and by multi- at ^ ^ u i- i. j ii, ■ -i 

•' ^ I Note that we have explicitly separated the cases j = j 

plying the final resuh by [l - Po{x)\ . This completes ^ ^ ^^^^^ take advantage of Eq. (|). We now 

the proof for K{x,x). again the decomposition ( p^ for /(x), [/(x)]^, and 

cr^ (the last taken to be a constant function of x\ thus 

2.1.2. Covariance obtaining 

Linearity has also important consequences for the form of ^ _ ^ ^ 

the covariance. In particular, the term (/(a;)/(a;')) can be (/(^)/(^ )) — / [/(a;)/(5;) -f (t ^C\[x^x \ x) 
written as 



(/»/>'))= / d5;p(5;)Ci(x,x';x)[a2 + /(x)/(x)] 



+ / Axp{x) \ Ax' p(x')j{x)!{x')C2{x,x';x,x') , (24) 



X 



with 

Ax p{x) j dx' p{x')C2{x,x']x,x')f{x)f{x' 



(19) Ci{x,x'-x) = -^^(^6{x - Xj)S{x- {{x,,5,j)}) 



Hence, the covariance is composed by a term proportional 
to cr^, the scatter of measurements {j/i} around the true x S{x'] {(i^i', <5i'j)}) ) , (25) 
values {/(a;^)} , and additional terms that can be identified / 

as Poisson noise. More specifically, the two noise terms are „ , / _ _/n 1 / ST^ c/- nc/-/ \ 

[cf.Eq.(§)] ^ C2{x,x-x,x) = -^^^{^5{x-x,)5{x - x,) 



T,{x,x') = J^dx p{x)Ci{x,x';x) , (20) x S{x; {ix„5,j)})S{x'; {{x^ ,5,,,,)})"^ 

Tp{x,x')^ / dx p{x)Ciix,x';x)[fix)]'^ 



(26) 



Again, the factors p{x) and p{x') have been introduced 
here to simplify some of the following equations. This 



r r nere to simpiiry some oi tne roiiowmg equations, ims 

+ j dxp{x) j dx' p(x')C2(x,x';x,x')f{x)f{x') proves Eq. (^I). Note that, although Eq. (H) is apparently 

composed of two independent factors, S{x\ ^{xi,5ij)^^ 



X Jx 



dx p{x) I dx' p{x')K{x;x)K{x';x')f{x)f{x') . and S(^x'; {{xii,Si>j>)]), in reality both terms are func- 



(21) 



tions of all locations {xi}; hence, since the random vari- 
ables that enters the expression for C2 are precisely the 



Again, the fact that measurement errors enter only locations {xi}, the two terms are correlated and Eq. (26) 
through CT^ is a consequence of linearity. cannot further simplified. The same, clearly, applies to 

In order to show Eq. ([T9|), we use again the linearity Eq. (p5|). 
of S to write [cf. Eq. ([l5|)] Similarly to the average, the two kernels Ci and C2 can 

be numerically evaluated using a simple procedure that we 
now describe: 



f{x) = S{x; {{xij{xi) + ei)}) 



N 



1. CI100S6 a, large subset A d X that contains x, 37, x\ 
^■^1 and x'. 
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2. Generate N, the number of objects inside A, according 
to the scheme chosen for the distribution of the loca- 
tions (for example, if the locations are independent and 
uniformly distributed with density p, then N follows 
the Poisson distribution given in Eq. (|^)). 

3. Generate the N locations {xi} inside A following the 
scheme chosen (for example, uniformly if the density 
is constant and there is no correlation). 

4. Assign a vanishing value yi = to each location Xi. For 
Ci, assign a value 1 to an extra location at x, similarly 
to what we do for K. For C2, assign a value 1 to two 
extra locations at x and x' . 

5. Hence, the two configurations considered are Ti = 
T = {(5,1)} U {(x„0)} and T2 = {(x, 1), (x', 1)} U 

6. For Ci evaluate S{x;Ti)S{x';Ti), if the function S 
is defined at both points x and x'; otherwise arbi- 
trarily assign a vanishing value to the product above. 
Analogously, for C2 evaluate S{x; T2)S{x'; T2), if this 
product is defined, or use zero otherwise. 

7. Repeat points 2-6, and evaluate the average 
(^(a;; Ti)S'(a;'; Ti)). This average, multiplied by [l — 

Po{x,x')] ^, is an estimate of Ci{x,x';x). An esti- 
mate of C2(x, x'; X, x') is instead given by the average 
{S{x- T2)S{x'- T2)) multiplied by [l - Po{x, x')] "\ 

Again, the numerical techniques described in the items 
above will be also used to derive some general properties 
for the covariance of interpolators in the following sections. 

This practical method to evaluate Ci and C2 can be 
derived using a technique similar to the one described for 
K. In particular, Eq. (|2^) suggests that we can determine 
the kernel Ci as follows. We approximate the delta distri- 
bution in X with a top- hat function \i[x — xj) of area s. 
Then, we randomly generate the locations inside a large 
subset A C X that contains x, x' , and x. As for the ker- 
nel K, most configurations will have only points outside 
the top-hat, thus resulting in a vanishing sum. A frac- 
tion sp{x) of configurations, however, will have a single 
point falling inside the top-hat function, i.e. with coor- 
dinate very close to x. Since the configurations without 
points inside s do not contribute to the average, we can 
just ignore them and consider only configurations with a 
point at x; we then multiply the average by the proba- 
bility of having such configurations, i.e. p{x)s in the limit 
s — > 0. Since the function H has a value 1/s at x, this is 
equivalent to assigning a value 1 to the point at x, and 
multiply the final result by p{x)] this term, however, dis- 
appears because of the presence of a factor l/p{x) in the 
definition of Ci. 

Again, in this procedure, we have ignored cases where 
the smoothing technique cannot be applied. In order to 
evaluate the expression in the right hand side of Eq. ( p5| ) , 
we need the values of 5 at x and x' . Hence, we should 
discard cases where S cannot be evaluated at x or x' . 
Alternatively, we can proceed as done for K, and just as- 
sign, for the configuration T, vanishing value to the prod- 
uct S{x;T)S{x';T) when one of the two factors cannot 



be evaluated. Finally, we multiply the total result by a 
factor [1 — Po{x,x')] , where this extra term is used to 
correctly normalize the average of Eq. (p5|). 

The evaluation for C2 is similar. In this case, how- 
ever, we have two different delta distributions at x and x' . 
Hence, we approximate them with two top-hat functions 
H(x — Xj) and H(a;' — x'j), both of area s. Since we are go- 
ing to take the limit s ^ 0, the two top-hat functions do 
not intersect, and thus the probability of having a point 
inside both top-hats is the product of the individual prob- 
abilities, i.e. p{x)p(x'). Then, in the limit s — > 0, we can 
force a point at x and one at x', both with value 1; other 
points are randomly distributed with density p. Then, we 
evaluate the average of 5'(x; T)5'(x; T) and multiply the 
final result by p{x)p{x'); this factor, however, disappears 
from C2 [cf. prefactor in Eq. (p^)] and is shifted instead 
in the integration (U). Again, if Pg is not vanishing, we 
need to correct for cases where S cannot be evaluated at 
X or x'; The correcting factor is still l/[l — Po{x, x')] . 

2.2. Normalization 

Almost every smoothing procedure used in Astronomy has 
a simple normalization property: If all values are the same, 
i.e. yi — c, then the interpolating function S always re- 
turns c at every point x: 

5(x;{(x„c)}) =c. (27) 

Normalized smoothing procedures are sometimes called 
"unbiased." This property is satisfied, for example, by the 
interpolator (^). For a linear smoothing, this property im- 
mediately implies 



K(x; x)p{x) dx = 1 



Vx e a: . 



(28) 



X 



This can be easily verified by using a constant function 
/(x) = 1 and by noting that in this case (/(x)) = 1; 
Eq. (|l|) then gives the normalization of K. Similarly, for 
the covariance kernels we have 

dx p(x)Ci(x, x'; x) 

X 

+ [ dxp(x) / dxV(x')C2(x,x';x,x') = 1 . (29) 
Jx Jx 

Note that the Poisson noise Tp vanishes if /(x) is "flat." 

Given any linear smoothing procedure S written in the 
form (^2]) , we can obtain a related normalized interpolator 
S': 



S"(x; {{x^,y^)}) = 



■ N 



N 



(30) 



X '^Wj{x;{x^}) 

Indeed, in many cases interpolating techniques are directly 
written as in Eq. (pO|). An example is given by our toy- 
interpolator (^, which is of the form ( pO| ) with 

Wj(x; {xi}) = l/jx - Xjl" . (31) 
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2.3. Spatial symmetries 

Basically all interpolation procedures of some interest sat- 
isfy some spatial symmetries. In this section we will con- 
sider three common spatial invariance properties, namely 
translation, rotation, and scaling invariance. These prop- 
erties can have important consequences on the forms of 
the three kernels K, Ci, and C2, provided that similar 
symmetries applies for the spatial distribution of location 
measurements . 

2.3.1. Invariance upon translation 

Many smoothing techniques are invariant upon transla- 
tion, i.e. there is no "preferred" point on X for the smooth- 
ing: 



S{x + d;{{xi + d,yi)}) = S{x; {{xi,yi)}) 



(32) 



For example, the simple smoothing (|2|) is invariant upon 
translation because it depends only on the distances 
\x — Xi\ between the interpolated point x and the loca- 
tions {xi}. A linear smoothing invariant upon translation 
has necessarily associated weights Wi [ci. Eq. (jl^)] of the 
form Wi = Wi[{xj — x}). Moreover, if the distribution 
of locations is also invariant upon translation (this im- 
plies, among other things, that the density p is uniform), 
then the kernel K is also invariant upon translation, i.e. 
K{x;x) = K{x — x), and thus Eq. (14) becomes a simple 
convolution. 

If the interpolation procedure is also normalized, than 
an interesting property holds: 



X 



{nx))dx^ / f{x)dx 



(33) 



X 



This can be shown by noting that for a smoothing proce- 
dure invariant upon translation, ^/^ is the convolution of 
/ with the kernel K. Moreover, for a normalized smooth- 
ing the kernel K is normalized to unity [cf. Eq. (p8|)], and 
thus Eq. (|3^ ) holds. Note that this result basically states 
the conservation of the "signal" : The total measured sig- 
nal is equal to the true, original signal. 

The invariance upon translation has also consequences 
on the two covariance kernels. More specifically, they 
must be of the form Ci{x,x';x) = Ci{x — x,x' — x) and 
C2{x, x'; X, x') — C2{x — x,x' — x; x! — x). Note also that 
Pq(x) is independent of x; hence, we will write just Pq; 
similarly, P{)(x,x') can depend only on [x — x'). 

2.3.2. Invariance upon rotation 

Smoothing procedures in spaces X of dimension larger 
than one are also often invariant upon rotation and mirror 
symmetry, i.e. they are isotropic. Formally 



S{x;{{x -I- R{xi - x),y^\) = S{x;{{xi,yi)\) 



(34) 



In this equation, i? is any orthogonal matrix; note that 
the rotation is actually performed around the point x. The 



smoothing (g) is invariant upon rotation (to show this we 
note, again, that this interpolator only depends on the 
distances \x — Xi\). 

If this symmetry holds also for the generation of lo- 
cations (this implies, again, that p is uniform) and if the 
interpolator is linear, than the kernel K is of the form 
K(x;x) = K' {x\ \x — xj), i.e. the smoothing kernel can 
depend on x and on the distance |a; — a;| between x and x 
only. 

The kernel associated with an interpolator invariant 
upon rotation has an interesting property. Let us evaluate 
the integral 



{x — x)K{x; x)dx — / (x — x)K(^x; \x — a;|) dx = , 
X Jx 

(35) 



where the last equality holds because the integrand is an 
odd function of (x — x). We can recast this result in a more 
interesting form 



xK(x; x) dx 



X 



K{x; x) dx 



X 



(36) 



This quantity in the left hand side of this equation is a 
measure of the "center of mass" for the kernel K; hence 



Eq. (36) basically assures that there is no systematic "off- 
set" on the final map. 

If both spatial invariance properties considered so far 
hold for a linear smoothing procedure, then the kernel 
K associated with the smoothing is of the form K{x; x) = 
K" {\x—x\) , so that only distances are involved. This prop- 
erty also holds for the covariance kernels, which are only 
dependent on the various distances between and 
(for C2 alone) x'. As a result, the noise term To- is a func- 
tion of |x — x'l alone. 

2.3.3. Invariance upon scaling 

A third spatial invariance occurs for some smoothing tech- 
niques, which are intrinsically scale-free: The result of the 
smoothing depends only on the ratios of distances and not 
on their absolute values. Formally in terms of the smooth- 
ing function, we have 



S{x;{{x + k{x, - x),yi)]) = S{x] {{x^,y^)]) , 



(37) 



where k is an arbitrary positive real number. It is not dif- 
ficult to verify that the interpolator (H) is invariant upon 
scaling. Indeed, for any positive real fc, we have 



S{x]{{x + k{x, - x),y,)]) 



E.L1/I- 



' Ef=ii/(fc"N-^.l") 

S{x-{{x,,y,)]) . (38) 



If measurement locations are uncorrelated and uni- 
formly distributed (i.e., for a homogeneous Poisson pro- 
cess), then the scale- free property allows us to derive a 
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number of consequences that greatly simplify the anal- 
ysis of the estimator. For a scale-free linear smoothing, 
the typical scale-length of the smoothing kernel K is set 
only by the density of objects p (no other scales are avail- 
able; note instead that in presence of clustering we would 
immediately have another scale, the correlation length). 
Formally, this can be seen by noting that two different lo- 
cation configurations, {xi\ and {a; -I- k[xi — x)} have the 
same probability if the density is changed from p to k~'^p, 
where n is the space dimension. As a result, all kernels are 
simply influenced by a change of p: 

p ^ k-''p , (39) 
K{x + d;x) ^ K{x + kd;x) , (40) 
Ci{x + d,x + d'; x) Ci{x + kd, x + kd'; x) , (41) 
C2{x + d, X + d' ; X , X + d) ^ 

C2{x + kd,x + kd';x,x + kd) . (42) 

Hence, we can study the kernels associated to a smooth- 
ing technique invariant upon scaling for a given density 
(say p = 1) and then generalize the results obtained using 
Eqs. (^-^^. We also note that for a normalized smooth- 
ing, as p cxD the product p{x)K will be proportional to 
Dirac's delta distribution. 

2.4. Exact interpolators 

An interpolating procedure is said to be exact if it honors 
the data points upon which the interpolation is based: 



S{xj; {{xi,yi))) = Vj 



(43) 



Honoring data points is seen as an important feature in 
many applications; on the other hand, if it is known that 
the data points might be affected by significant errors, 
exact interpolators might not be the right choice. 

Regarding the exactness, the interpolator ^ needs a 
few words of explanation. Indeed, when the interpolated 
position X coincides with one of the locations {xi\, the 
expression (|^) becomes singular, and thus the interpolator 
is apparently undefined. However, the limit x — > of S" is 
defined and equal to yi (we assume here that all locations 
are different; the case of two or more coincident locations 
has no statistical relevance). Hence, we can safely extend 
the definition of the smoothing technique (||) to all points 
x & X, with the convention that limits must be taken if 
the expression is undefined. This way, we can state that 
our toy- interpolator is exact. 

Recalling the numerical method used to evaluate K 
for a linear interpolator, we can deduce that an exact es- 
timator has the property K{x; x) = 1/ [l — ^o(x)] • In fact, 
if we use the numerical method, we will always measure 
S'(x; {(xi, ?/i)}) = 1, since we have assigned a value 1 to 
the location x, and thus the average of all those measure- 
ments is also 1 (here we are neglecting cases where the 
smoothing is not defined). Similarly, we find Ci[x,x;x) — 
l/[l — Po(x,x)] and C2{x,x';x,x') — C2{x' , x; x , x') = 
l/[l-Po{x,x')]. 



2.5. Bounded interpolators 

An interpolating procedure is said to be bounded if inter- 
polated values are always between the smallest and the 
largest measured values: 



iyi < S{x; {{xi,yi)}) < maxy^ 



(44) 



This is a rather natural property satisfied by many 
smoothing techniques [e.g., the interpolator (||)]. Note 
that a bounded interpolator is always normalized, since 
Eq. (H) implies Eq. @). 

If a linear interpolator is bounded, then we can put 
superior and inferior limits on K{x;x). In fact, recalling 
again the numerical technique used to evaluate K , we ob- 
tain 



< K{x;x) < [l-Po{x)\ 



(45) 



Hence, the kernel K is non-negative and [l — ^o(x)] ^ is 
a superior limit for it. Similarly, for the other kernels we 
obtain 



0<Ci{x,x';x)<[l-Po{x,x')] \ 
< C2{x, x ; X, x) < [l — Pq{x, x')\ ^ . 



(46) 
(47) 



We can go further in this analysis by defining the effec- 
tive number of objects used by the smoothing procedure 
as 



Ks{x) = 



1 



l-Po(x) 



K{x; x)p{x) dx 



X 



[K{x;x)\ p{x)Ax 



X 



(48) 



Note that A/'eft can be defined for any linear smooth- 
ing procedure; for a normalized (and in particular for a 
bounded) interpolator, the first term in the right hand 
side of Eq. (^8|) is unity. For an interpolator invariant 
upon translation, the quantity 7Vcfr(x) does not depend 
on X, and thus we will write just TVog. Using the defini- 
tion (^) we can recast the upper limit for K{x;x), valid 
for a bounded interpolator, as a lower limit for the effec- 
tive number of objects A/'cff (x): 



A4ff'(x) > 



1 



1 

> 1 , 



Po(x) 



1 



X 



l-Po(x 



■K{x; x)p(x) dx 



(49) 



Equation ( p^ ) put a lower limit to the effective number for 
the kernel K , i.e. to the expected resolution of the smooth- 
ing. Finally, we observe that interpolators for which the 
scaling invariance holds have effective number indepen- 
dent of p (remember that p is supposed to be constant for 
scaling invariant interpolators). 

The covariance of a bounded interpolating procedure 
satisfies some interesting properties. Using Ci > [see 
(|6|)] in Eq. (pi), in fact, we obtain 



dx p(x)Ci(x, x'; x) < 1 , 



(50) 



X 
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so that Trj{x,x') < a^. Hence, we have obtained an upper 
limit for the contribution to the noise from measurement 
errors. A comparison of this result with the inequahty 
■Miff > 1 suggests that our definition of effective number 
of objects is sensible. 

We can better appreciate the relationship between the 
covariance kernel Ci and MeS with the following argu- 
ment, valid for any linear interpolation procedure. Let 
us consider Ci{x,x]x) and briefly recall how this quan- 
tity can be measured using the numerical approach dis- 
cussed in Sect. 2.1.2 . We generate locations, keeping one 
object with weight 1 in x, and take the average of the 
quantity 5'(x; Ti)S'(a;; Ti); finally we multiply the aver- 
age by l/[l — Po(2;)], because in this particular case 
Po{x,x) = Pq{x). Something quite similar is done for 
K{x]x) (see Sect. ^.1.1| ): The configurations used are the 
same, in the sense that T = Ti, but in this case we take 
the average of S{x; T). From this discussion, we see that 
K(x; x) is the simple average of a set of some non-negative 
numbers, while Ci{x, x; x) is, apart from a known numer- 
ical factor, the average of the squares of the same quanti- 
ties. Hence, the following inequality holds: 



Ci{x,x;x) > [l-Po{x)] [K{x;x)y 



(51) 



For a normalized interpolator, this equation implies 
Ta-{x, x) > (T^/A/"cff . Hence, the effective number of objects 
puts a lower limit on the measurement error. 

2.6. Local interpolators 

An interpolator is called local if the interpolated value at 
x depends only on the values yi for locations Xi "close" to 
x; an interpolator which is not local is said to be global. 
Global interpolators usually produce smoother maps but 
often require more computing time. The interpolator (^ 
is global. 

We further distinguish between strongly and weakly lo- 
cal interpolators depending on the meaning that is given 
to the word "close" in the previous definition. An inter- 
polator is strongly local if it uses only values jji whose 
locations Xi fall within a fixed range from x; a weakly lo- 
cal interpolator, instead, is only guaranteed to use a finite 
number of points close to x. We stress that this definition 
is not standard but is needed for a rigorous characteriza- 
tion of local smoothing techniques. Note that a strongly 
local interpolator cannot be scale invariant, because it de- 
pends only on points within a fixed range and hence, there 
is a natural scale for the interpolator, its range. A strongly 
local interpolator has always Po{x) > 0, because the inter- 
polator cannot be defined if there is a large "void" without 
points around x. 

Using again the numerical technique to obtain K, 
we immediately find that the kernel K associated with 
a strongly local interpolator has compact support, i.e. 
K{x; x) vanishes if ja: — 5;| is large; analogously, Ci(x, x'; x) 
vanishes if either \x — x\ or \x' — x\ are large, and 
C2(x, x'; X, x') vanishes if x or x' are far away from x or 



x'. Note also that for strong local estimators, if \x — x'\ is 
large, C2{x,x';x,x') converges to 



C2{x, x'; X, x') 



[l-Pn[x)] [l-P,{x')\ 



[l~Po{x,x')\ 
X [K{x]x)K{x']x')+K{x';x)K{x]x')\ . (52) 

Not much, instead, can be said regarding weakly local in- 
terpolators; indeed, weakly local interpolators have been 
defined here mainly to have a complete classification of 
locality. 

Several strongly local interpolators can be defined with 
just a single point inside their range. In this case, we can 
show that in the limit p ^ 0, the kernel K associated to a 
interpolator approaches a top-hat function, provided the 
density p is uniform. For if p is sufficiently small every- 
where, we expect a negligible probabihty, when using the 
numerical technique to obtain K{x]x)^ of having points 
other than x inside the range of the function. In this case, 
all terms in the sums of Eq. (30) vanish except for x. As 
a result, we always measure 1 for points inside the range 
of the interpolator, and outside. Moreover, if tt is set 
around the point x where the interpolator is defined, we 
find [cf. Eq. (H)] 

Po{x) = cxp^- j p{x') dx'^ - J • 

Notice that the dependence on x of K{x; x) is also ex- 
pressed by the set tt, which is a function of the location x. 
In summary, we obtain 



K{x; x) 



1/ p{x') Ax' ff X G TT , 



otherwise 



(54) 



As expected, this kernel is manifestly normalized accord- 



ing to Eq. (27). In case of a uniform distribution p, this 
expression simplifies to 



K{x; x) 



1/ PP{it) if a; G tt , 
otherwise , 



(55) 



where p{'k) is the area of the set tt. Hence, in the limit 
p ^ we reach the upper bound for K stated in Eq. ( ^ . 
Regarding the covariance kernel, calling /i(7rn7r') the area 
of the intersection of the ranges for x and x' , we obtain 
(here, for simplicity, we assume a uniform location density 
P) 



Ci{x, x' ; x) 



I l/p/x(7r n tt) if a; e tt n tt 



1° 

6*2(2;, x'; X, x') ~ 



otherwise 



(56) 
(57) 



Note that C2 vanishes in the limit p ^ (this is due to the 
extra factor p, with respect to Ci , present in the evaluation 
of C2). Note also that in this case we have reached the limit 
Ta — discussed after Eq. (pO|). Finally, we can verify 
without difficulties that the normalization of Eq. ( p9[ ) is 
satisfied. 
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2. 7. Other properties 

Recalling again the numerical technique to obtain K de- 
scribed in Sect. |2.1.1| , we can consider a number of other 
classifications all focused on properties of S when all but 
one point are vanishing. 

For example, suppose that the smoothed map S for 
the configuration T = {(x, 1)} U {{xi,0)} is smaller than 
unity at each point. Such a smoothing procedure could 
not be bounded, but we can still obtain the upper limit 
K{x;x) < 1/ (1 — Pq{x)^ ioT K{x;x) (as in case of bounded 
interpolators) and Eq. (^) holds. 

Another interesting property to consider is monotonic- 
ity. Suppose that, for the same configuration T consid- 
ered above, the function S is monotonically decreasing as 
\x — x\ increases: Then K(x;x) will also decrease mono- 
tonically with |a; — xj. If the interpolator is exact, then, 
this property implies K{x;x) < l/[l — PqC^^)]; because 
K{x, x) — 1/ [l — ^b(a;)] , and again we recover Eq. (^. 

3. Kernels of some smoothing procedures 

In this section we will consider several smoothing proce- 
dures and derive, analytically or numerically, the relative 
kernels K; for simplicity, we will not consider the kernels 
Ci and C2; moreover, we will assume a uniform spatial 
distribution of locations characterized by a density p and 
no correlation. A standard approach will be followed for 
each interpolator. We will first briefly introduce and define 
the interpolator; then we will classify it according to the 
nomenclature introduced in Sect. ^ finally we will evalu- 
ate the kernel K. 

The interpolation techniques discussed below are gen- 
erally well known and are described in any book on spa- 



tial interpolation (e.g. Watson 1992). As a reference (ex- 
tremely useful also for a deep discussion on statistics of 



spatial data) we refer to Cressie (1993), and in particu- 



lar to Sect. 5.9 of this book (see also Adlcr 1991; Preston 
197^ ). 



3.1. Nearest neighbor 

This interpolator, called also "proximal," is probably the 
simplest one can imagine: The value in a point is assumed 
to be equal to the value of the nearest known point. Hence, 
the output data for such interpolator consists in Voronoi 
cells (see Okabe et al. 1992) centered on the points {xi} 



with abrupt changes at the boundaries (see Fig. 0). 

This interpolator is manifestly linear, normalized, 
invariant upon translation, rotation, and scale, exact, 
bounded, and weakly local. As a result, we expect a ker- 
nel K of the form K{x;x) = K (\x — x\) which scales with 
the density according to Eqs. ( ^ and (pO|). Moreover, 
we expect Pq{x) = Pq(x,x') — 0, AfcS independent of p, 
K{0) = 1, and K{r) < 1. Finally, K{r) will be not in- 
creasing. 

For this estimator, we can actually carry out cal- 
culations analytically using the approach discussed in 




Fig. 1. The relationship between Voronoi cells (solid lines) 
and Delaunay triangles (dashed lines). The figure has 
been made by calculating the Delaunay triangulation and 
Voronoi tessellation for a large set of points on the plane, 
and by showing a small region. Note that two points are 
connected in the Delaunay triangulation if and only if their 
Voronoi cells have a side in common. 



Sect. |T|. Let us consider a point x at distance r from 
X. The probability that the closest location for x be a; is 
given by p{r) — exp(— pfc„r"), where n is the dimension 
of X and kn is the measure (area or volume) of the unit 
sphere (e.g., ki = 2, k2 — tt, and fcs — Att/3). In fact, x 
is the closest location to x if the ball of radius r centered 
on X does not contain any other location. The number 
of locations inside this ball follows a Poisson distribution 
with average pfc„r", and thus this number vanishes with 
the probability p{r) given above. Hence, the point at x 
will have a value 1 with probability p{r), and a vanishing 
value with probability 1 — p(r). In summary we obtain 



K(^\x — x\) — exp(— pfc„|a; — x\") 



(58) 



Using this expression for K we can check all the properties 
stated in Sect. I (cf. also Figs. I and |). Finally, we have 
A4ff = 2. 



3.2. Delaunay triangulation 

The nearest neighbor technique is very simple but unfor- 
tunately produces abrupt changes at the boundaries of 
the Voronoi cells, i.e. this interpolator is not smooth. This 
problem can be solved by using n+l near points (where n 
is the space dimension) and by linearly interpolating their 
values. For example, if X = R, so that n = 1, we take the 
point at the left and the point at the right of x and then 
use linear interpolation. If n > 1, however, we have the 
problem of choosing the most convenient points to perform 
the linear interpolation (mathematically speaking, R" is 
not totally ordered for n > 1). A standard approach is 
to use Delaunay triangulation (see, e.g. p' Rouche 1994), 
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Fig. 2. The kernels K associated with the one-dimensional 
nearest neighbor (solid line) and Delaunay interpolators 
(dashed line) for p = 1. The kernel K for the natu- 
ral neighbor smoothing is identical to the one for the 
Delaunay interpolator. Since the interpolators considered 
are invariant upon translation and scaling, the kernels as- 
sociated to different densities can be evaluated using the 
transformation ( ^9|) and (^0|). Note also that for both plots 
we find K{0) — 1 and K{x; x) <1 because the smoothing 
techniques are exact and bounded. 



that we describe here for n = 2. A Delaunay triangula- 
tion for the locations {xi} is the unique triangulation of 
the plane with the property that no point in the set {xi} 
falls in the interior of the circumcircle (circle that passes 
through all three vertices) of any triangle in the triangu- 
lation. This definition is easily generalized for n> 2. The 
Delaunay triangulation for the points {xi\ is also closely 
related to the Voronoi cells for the same points, one being 
the dual of the other. For example, the Delaunay trian- 
gulation can be obtained by connecting all points whose 
Voronoi cells have a side in common (see Fig. |l|). 

The Delaunay triangulation interpolator is thus ob- 
tained in the following way. The Delaunay triangulation 
is calculated for the points {xi}] the interpolated value at 
X is then obtained using linear interpolation of the values 
for the three points that are the vertices of the triangle to 
which X belongs. Delaunay interpolation, or variants of it, 
has found already many applications in the astrophysical 



context (see, e.g. 


van de Weygaert 


1994 




Bernardeau & 


van de 


Weygaert 


1996; 


Schaap & van de Weygaert 


2000|). 



The Delaunay triangulation interpolator is linear, nor- 
malized, invariant upon translations, rotations, and scal- 
ing, exact, bounded, and weakly local. We expect for K 
the same properties as for the nearest neighbor. 

Analytical calculations for the Delaunay triangulation 
linear smoothing are trivial only for dimension n — 1. In 
this case, in fact, we can just use the point at the left 
and the point at the right of x and proceed with linear 
interpolation. Calling d_ > and d+ > the distances 
between the closest locations at the left and at the right 
of X and x itself, we have 



d_f{x + d+) + d+f{x - d_) 



The probability distribution for both d+ and d- is 
p{d±) = pe-P''^ , 
and thus we obtain 



(60) 



{f{x)) = P~ I C • ' UU+ 
'0 



oo 



X / e 
/o 

/•oo 



p-^+dd- 

_^^_f{x + d+)d^+f{x^d^)d+ 



d+ + rf_ 
p2 / f{x + xi)e-P^''^Uxi 



dd_ 
dx2 



\xi\+X2 

p I f{x)K{x ~ x)dx , (61) 



where the kernel K is given by 

Kix) = e-^'l^l - p|a;|r(0, p\x\) . (62) 
Here F is the incomplete gamma function, defined as 



r(a, x) 



(63) 



The sa me res ult is obtained using the numerical technique 
of Sect. |2.1.l| . It can be verified that, as expected, K{0) = 
1 and K{x) < 1; moreover, K depends only on the product 
pI^I, so that the scaling invariance (|39|-|4C|) holds. For this 
kernel we find Wcfr ~ 2.44417. A plot of K{x) for n = 1 is 
shown in Fig. ^. 

For 71 > 1 calculations cannot be easily carried out 
analytically. However, since the smoothing procedure is 
expected to be scale invariant, we can numerically evalu- 
ate K for, say, p — I, and then apply Eqs. ( ^ and (40) 
to convert this result to different densities. Note also that, 
because of the translation and rotation invariance proper- 
ties, K is expected to be of the form K(\x — ij), i.e. it is 
a function of a single real argument. Figure |^ shows the 
numerical results obtained for dimension n = 2; again, we 
can check that all expected properties for K are indeed 
satisfied. For the two-dimensional Delaunay interpolator 
we find AAeff — 3.31. 

3.3. Natural neighbor 



(59) 



The natural neighbor interpolator (Sibson 1981 ), also 
called area stealing interpolator, is still related to the 
Voronoi cells. Its construction is carried out in the fol- 
lowing way. First, the Voronoi cells for the points {xi} are 
calculated. Then, the point x is added to the set, and the 
new Voronoi cells are calculated. The cell corresponding 
to X in the new configuration overlaps some parts (stolen 
areas) of cells originally owned by nearby points: These 
points, called natural neighbors, will be involved into the 
interpolation at x. Specifically, the value assigned to x 
will be the weighted average [see Eq. (|3^)] of values at the 
natural neighbors, with weights equal to the overlap areas 
(see Fig. |). 

The natural neighbor interpolator is linear, normal- 
ized, invariant upon translation, rotation, and scaling, 
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Fig. 3. The kernels K associated with the two- 
dimensional nearest neighbor (solid line), Delaunay trian- 
gulation (dashed line) , and natural neighbor (dotted line) 
interpolators. For all plots we have used a density p ~ 1. 
Since all smoothing techniques are invariant upon trans- 
lation, rotation, and scaling, the kernels associated to dif- 
ferent densities can be evaluated using the transformation 
(H) and dH). 



area will almost be entirely inside the (old) cell corre- 
sponding to X3. In the limit where x coincides with 2:3, 
the stolen area will be the intersection of the cell of X3 
with the half-plane {x \ (x — X3,x) > O}, where (•, •) 
denotes the scalar product. 

Again, analytical calculations to obtain K for this 
smoothing technique are feasible only in dimension n — I. 
In this case, surprisingly, the natural neighbor smoothing 
is totally equivalent to the Delaunay triangulation, and 
thus Eq. (||) holds. 

For higher dimensions, wc can only obtain numerical 
estimates for iir(|x — x|) . Figure || shows these estimates 
in dimension n = 2; in this case we obtain AfcS — 3.41. 

3.4. Moving weights 

One of the most common interpolating techniques used in 
Astronomy is based on a simplified form of Eq. (^^ in 
which the i-th weight depends only on x and on the i-th 
location: 




r N 



S{x; {{x^,y^)}) = 



r N 



'^w{x-Xi)yi / y^^wjx - Xj) 



Fig. 4. The evaluation of coefhcients for the natural neigh- 
bor interpolation. A first Voronoi tessellation for the orig- 
inal set of locations (empty circles) is calculated. Then, 
the point x (filled circle) for which we want to obtain the 
interpolated value is added to the original locations, and 
a second Voronoi tessellation is evaluated. The new tessel- 
lation coincides with the old one, with the only difference 
that the new point x is now "stealing" the area for its 
Voronoi cell from neighboring cells (the ones correspond- 
ing to the points Xi, X2, x^, and X4 in the case considered 
in this figure). The various "stolen areas" (marked in this 
figure with parallel lines of different orientation) are then 
used to evaluate the coefficients to be inserted in Eq. (BQ). 



exact, bounded, and weakly local. All these properties 
are easily verified except, perhaps, the exact property. 
Suppose that the point x moves toward one of the loca- 
tions, say ^3. As a; get closer and closer to 0:3, the stolen 



(64) 



The function wlx — Xi), taken here to be non-negative, is 
usually chosen to have a peak aX x = Xi and to decrease to 
zero as |a; — x^j increases. This way, the estimated value at 
X will be basically determined by the neighboring points. 
Often the weight function is isotropic, i.e. w{x) — w{^\x^. 
Commonly used weight functions include Gaussians, top- 
hats, and inverse distances of the form 'w{x) = l/|a;|" 
with a fixed [this is precisely our toy-interpolator defined 
in Eq. (1)]. 

Some of the properties of the moving weight interpola- 
tor strongly depend on the weight function used. In gen- 
eral, this interpolator is linear, normalized, invariant upon 
translations, and bounded. Hence, the kernel associated 
to the smoothing satisfies K{x — a;) < 1 and is normal- 
ized to unit. The interpolator is invariant upon rotation 
if and only if w is isotropic. If w is homogeneous of de- 
gree g (as for example for inverse distances weights), i.e. if 
w{kx) — k^w{x), then the interpolator is invariant upon 
scaling . The interpolator is generally not exact; however, 
if w{x) is bounded for x ^ Q and goes to infinity for a; = 
(e.g., in case of inverse distance weights), then the inter- 
polator is exact and thus we expect K{0) = 1. The inter- 
polator is strongly local if and only if wlx) has compact 
support, in which case K has also compact support; it 
is otherwise global. Finally, we observe that this interpo- 
lator is defined if just one point is inside the support of 
w. Hence, if w has compact support of area a, we have 
Po{x) = exp[— pa]; in other cases Po{x) vanishes. 

The statistical properties of this smoothing have been 
studied in detain in a separate paper (Lombardi & 
iSchncidcij ^001\ ). There we have shown that the kernel 
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can be evaluated using the set of equations 
K{x) 



-sw{x) 



1 



X 

w{x) 



dx 



1-Pr 



-w{x)s+pQ{s) 



Jo 



(65) 
(66) 



In the same paper we have also explicitly shown several of 
properties stated above for this interpolator. In addition, 
we have shown that the kernel K associated with this 
smoothing closely resemble the weight function w used, 
but is generally larger. For example, if we define the area 
Af of w analogously to Eq. (||), we find AfeS > N. The 
kernel K however converges to w as p — > cxo. 

3.5. Fixed weights 

For some studies a non-normalized form of the smoothing 
(Q) is still used: 



can be found in de Boor 1978). They find several applica- 
tions in many different fields in Astronomy, and can also 
be used for "cosmetic treatments." Splines come in sev- 
eral variants with slightly different properties. For exam- 
ple, some variants of splines are exact interpolators, while 
others are not. Another peculiarity is that often splines 
are not bounded interpolators. In order to avoid any pos- 
sible confusion among the different variants of splines, we 
prefer here not to consider this smoothing technique in 
detail. We stress again, however, that a statistical anal- 
ysis of any particular variant can be carried out without 
difficulties using the framework described in this paper. 

Kriging is another interesting interpolation method. 
Although it is seldom used in Astronomy, it finds sev- 



eral applications in geophysics (see, however, [Alfaro et al 



JV 

S{x; {{xi.yi)]) = ^w{x - Xi)y, 

i=l 



(67) 



1991 for an example of astronomical application of krig- 
ing). A complete discussion of this interpolation method, 
which also has some variants, is beyond the scope of this 
paper. Here, we just observe that kriging is a non-linear 
method, and thus cannot profitably be analyzed with the 
techniques discussed in this article. 



This smoothing has been employed, e.g., in some weak ^- Conclusions 



lensing studies ( 


Kaiser & Squires 


1993; 


Tyson & Fischer 


1995; Luppino & Kaiser 


1997; 


Fischer 


199C). 



properties described above. In particular, it is only lin- 
ear and invariant upon translation. If w(x) = w{^x{) is 
isotropic, then the rotation invariance holds. The inter- 
polator is not normalized (which is quite unusual for a 
smoothing technique); however, if the weight function w 
is normalized to 1, i.e. 



In this paper we have considered the statistical properties 
of interpolation techniques. The discussion has originally 
been kept at a high level of abstraction, and this has given 
us the ability to characterize smoothing methods in terms 
of simple properties. In particular, we have made a clas- 
sification of interpolation techniques and we have derived 
several statistical results associated with each class of in- 
terpolators. A comparison of this analysis with a more 
technical one carried out in a separate paper for a specific 



p I w{x) dx 
'x 



1 



(68) 



then Eq. (g8|) still holds. The interpolator is not exact 
and not bounded; it is strongly local if w has compact 
support, otherwise it is global; it is always defined, so that 
Poix)^PQ{x,x')^0. 

A statistical analysis of this interpolator is straightfor- 
ward. The result obtained is that the kernel K coincides 
with the weight function w. As expected, the kernel is not 
necessary bounded su periorly and is not ne cessary equal 
to unit at x = (see [Lombardi et al. 2001 for a detailed 
discussion of this smoothing technique applied to weak 
lensing data). 

3.6. Other interpolation techniques 

Clearly, a complete analysis of all smoothing methods cur- 
rently used would be impossible here; on the other hand, 
the general classification technique described in Sect. ^ 
should allow one to characterize other interpolators. We 
want to mention however a couple of interesting smooth- 
ing techniques, namely the splines and the kriging. 



Splines are an interesting and complex topic (see Press 
et al. 1992| for a short introduction; a complete discussion 



smoothing method (Lombardi & Schneider 2001) clearly 
shows the advantages of using an abstract approach to the 
problem. In the second part of this paper we have applied 
the results obtained to several commonly used smoothing 
methods. 
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